---
layout: post
title: "Log-likelihood benchmark"
date: "2015-09-08"
author: Alex Rogozhnikov
tags:
- notebook
---


<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Log-likelihood-benchmark">Log-likelihood benchmark<a class="anchor-link" href="#Log-likelihood-benchmark">&#182;</a></h2><p>This is a simple benchmark, which I use for basic test of vector-based computing engines.</p>
<p>The problem is to find log-likelihood of normal distribution:
$\sum_i \log \left[ \dfrac{1}{\sqrt{2 \pi \sigma}} \exp \left( - \dfrac{(x_i - x)^2}{2 \sigma^2} \right) \right] $</p>
<p>There are elementwise subtraction, division, power, exponent, logarithm and summation of array, so this is broader test.</p>
<p>tested:</p>
<ol>
<li>numpy + scipy's pdf</li>
<li>numpy</li>
<li>cython</li>
<li>numexpr</li>
<li>theano</li>
<li>parakeet</li>
<li>C++</li>
<li>FORTRAN</li>
</ol>
<p>computing log-likelihood for normal distribution</p>
<p><strong>Notes</strong></p>
<ol>
<li>Not optimizing computations here (but in theory theano and parakeet may remove unnecessary computations)</li>
<li>This test includes exp, log, division and summation of array</li>
<li>Everything is running on CPU in one thread, this limitation is for 'fairness' of tests</li>
</ol>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[1]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">numpy</span>
<span class="kn">import</span> <span class="nn">scipy</span>
<span class="kn">from</span> <span class="nn">scipy.stats</span> <span class="kn">import</span> <span class="n">norm</span>
<span class="kn">import</span> <span class="nn">theano</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Scipy-implementation">Scipy implementation<a class="anchor-link" href="#Scipy-implementation">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[2]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="k">def</span> <span class="nf">llh_scipy</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="n">lh</span> <span class="o">=</span> <span class="n">norm</span><span class="p">(</span><span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">)</span><span class="o">.</span><span class="n">pdf</span><span class="p">(</span><span class="n">data</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">numpy</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">lh</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Numpy-implementation">Numpy implementation<a class="anchor-link" href="#Numpy-implementation">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[3]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="k">def</span> <span class="nf">llh_numpy</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="n">s</span> <span class="o">=</span> <span class="p">(</span><span class="n">data</span> <span class="o">-</span> <span class="n">mean</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span> <span class="o">/</span> <span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="p">(</span><span class="n">sigma</span> <span class="o">**</span> <span class="mi">2</span><span class="p">))</span>
    <span class="n">pdfs</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span> <span class="n">s</span><span class="p">)</span>
    <span class="n">pdfs</span> <span class="o">/=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">pi</span><span class="p">)</span> <span class="o">*</span> <span class="n">sigma</span>
    <span class="k">return</span> <span class="n">numpy</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">pdfs</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Cython-implementation">Cython implementation<a class="anchor-link" href="#Cython-implementation">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[4]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">%</span><span class="k">load_ext</span> Cython
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[5]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">%%</span><span class="n">cython</span>
<span class="k">cdef</span> <span class="kr">extern</span> <span class="k">from</span> <span class="s">&quot;math.h&quot;</span><span class="p">:</span>
    <span class="n">double</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">double</span> <span class="n">m</span><span class="p">)</span>
    <span class="n">double</span> <span class="n">exp</span><span class="p">(</span><span class="n">double</span> <span class="n">m</span><span class="p">)</span>
    <span class="n">double</span> <span class="n">log</span><span class="p">(</span><span class="n">double</span> <span class="n">m</span><span class="p">)</span>

<span class="k">import</span> <span class="nn">cython</span>
<span class="k">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>

<span class="k">cimport</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="k">from</span> <span class="nn">numpy</span> <span class="k">cimport</span> <span class="n">ndarray</span>
<span class="n">pi</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span>

<span class="nd">@cython</span><span class="o">.</span><span class="n">boundscheck</span><span class="p">(</span><span class="bp">False</span><span class="p">)</span>
<span class="nd">@cython</span><span class="o">.</span><span class="n">wraparound</span><span class="p">(</span><span class="bp">False</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">llh_cython</span><span class="p">(</span><span class="n">ndarray</span><span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">float64_t</span><span class="p">,</span> <span class="n">ndim</span><span class="o">=</span><span class="mf">1</span><span class="p">]</span> <span class="n">data</span><span class="p">,</span> <span class="n">double</span> <span class="n">mean</span><span class="p">,</span> <span class="n">double</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="k">cdef</span> <span class="kt">int</span> <span class="nf">l</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">data</span><span class="p">)</span>
    <span class="k">cdef</span> <span class="kt">double</span> <span class="nf">llh</span> <span class="o">=</span> <span class="mf">0</span>
    <span class="k">cdef</span> <span class="kt">double</span> <span class="nf">s</span> <span class="o">=</span> <span class="mf">0</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">l</span><span class="p">):</span>
        <span class="n">s</span> <span class="o">=</span> <span class="p">(</span><span class="n">data</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span> <span class="n">mean</span><span class="p">)</span> <span class="o">**</span> <span class="mf">2</span> <span class="o">/</span> <span class="p">(</span><span class="mf">2</span> <span class="o">*</span> <span class="p">(</span><span class="n">sigma</span> <span class="o">**</span> <span class="mf">2</span><span class="p">))</span>
        <span class="n">llh</span> <span class="o">+=</span> <span class="n">log</span><span class="p">(</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">s</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2</span> <span class="o">*</span> <span class="n">pi</span><span class="p">)</span> <span class="o">*</span> <span class="n">sigma</span><span class="p">))</span>
    <span class="k">return</span> <span class="n">llh</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Numexpr-implementation">Numexpr implementation<a class="anchor-link" href="#Numexpr-implementation">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[6]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">numexpr</span>
<span class="n">numexpr</span><span class="o">.</span><span class="n">set_num_threads</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">llh_numexpr</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="n">expression</span> <span class="o">=</span> <span class="s">&#39;sum(log(exp(- (data-mean) **2  / (2 * sigma ** 2)) / (sqrt(2 * pi) * sigma)))&#39;</span>
    <span class="k">return</span> <span class="n">numexpr</span><span class="o">.</span><span class="n">evaluate</span><span class="p">(</span><span class="n">expression</span><span class="p">,</span> <span class="n">local_dict</span><span class="o">=</span><span class="nb">dict</span><span class="p">(</span><span class="n">data</span><span class="o">=</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="o">=</span><span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="o">=</span><span class="n">sigma</span><span class="p">,</span> <span class="n">pi</span><span class="o">=</span><span class="n">numpy</span><span class="o">.</span><span class="n">pi</span><span class="p">))</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Parakeet-implementation">Parakeet implementation<a class="anchor-link" href="#Parakeet-implementation">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[7]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">parakeet</span>
<span class="n">parakeet</span><span class="o">.</span><span class="n">config</span><span class="o">.</span><span class="n">backend</span> <span class="o">=</span> <span class="s">&#39;c&#39;</span>

<span class="nd">@parakeet.jit</span> 
<span class="k">def</span> <span class="nf">llh_parakeet</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="n">s</span> <span class="o">=</span> <span class="p">(</span><span class="n">data</span> <span class="o">-</span> <span class="n">mean</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span> <span class="o">/</span> <span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="p">(</span><span class="n">sigma</span> <span class="o">**</span> <span class="mi">2</span><span class="p">))</span>
    <span class="n">pdfs</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span> <span class="n">s</span><span class="p">)</span>
    <span class="n">pdfs</span> <span class="o">/=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">pi</span><span class="p">)</span> <span class="o">*</span> <span class="n">sigma</span>
    <span class="k">return</span> <span class="n">numpy</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">pdfs</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Theano-implementation">Theano implementation<a class="anchor-link" href="#Theano-implementation">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[8]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">theano</span>
<span class="kn">import</span> <span class="nn">theano.tensor</span> <span class="kn">as</span> <span class="nn">T</span>

<span class="k">print</span> <span class="n">theano</span><span class="o">.</span><span class="n">config</span><span class="o">.</span><span class="n">device</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>cpu
</pre>
</div>
</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[9]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">theano</span><span class="o">.</span><span class="n">config</span><span class="o">.</span><span class="n">openmp</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[9]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>False</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[10]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="k">def</span> <span class="nf">llh_theano</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="n">s</span> <span class="o">=</span> <span class="p">(</span><span class="n">data</span> <span class="o">-</span> <span class="n">mean</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span> <span class="o">/</span> <span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="p">(</span><span class="n">sigma</span> <span class="o">**</span> <span class="mi">2</span><span class="p">))</span>
    <span class="n">pdfs</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span> <span class="n">s</span><span class="p">)</span>
    <span class="n">pdfs</span> <span class="o">/=</span> <span class="n">T</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">pi</span><span class="p">)</span> <span class="o">*</span> <span class="n">sigma</span>
    <span class="k">return</span> <span class="n">T</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">pdfs</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>

<span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">scalars</span><span class="p">(</span><span class="s">&#39;m&#39;</span><span class="p">,</span> <span class="s">&#39;s&#39;</span><span class="p">)</span>
<span class="n">d</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">vector</span><span class="p">(</span><span class="s">&#39;data&#39;</span><span class="p">)</span>

<span class="n">llh_theano</span> <span class="o">=</span> <span class="n">theano</span><span class="o">.</span><span class="n">function</span><span class="p">([</span><span class="n">d</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">],</span> <span class="n">llh_theano</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">))</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="FORTRAN-implementation">FORTRAN implementation<a class="anchor-link" href="#FORTRAN-implementation">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[11]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">%</span><span class="k">load_ext</span> fortranmagic
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>


<div class="output_subarea output_javascript ">
<!--<script type="text/javascript">
$.getScript("https://raw.github.com/marijnh/CodeMirror/master/mode/fortran/fortran.js", function () {
IPython.config.cell_magic_highlight['magic_fortran'] = {'reg':[/^%%fortran/]};});

</script>-->
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[12]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">%%</span><span class="k">fortran</span> 
subroutine llh_fortran(data, mean, sigma, result)
    real*8, dimension(:), intent(in) :: data
    real*8, intent(in) :: mean, sigma
    real*8, intent(out) :: result
            
    real*8, dimension(size(data, 1)) :: s
    real*8, parameter :: PI = 3.14159265358979323846

    s = (data - mean) ** 2 / (2 * sigma ** 2)
    s = exp(- s) / (sqrt(2 * PI) * sigma)
    result = sum(log(s))
    
end subroutine llh_fortran
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="C++-implementation-for-comparison-of-speed">C++ implementation for comparison of speed<a class="anchor-link" href="#C++-implementation-for-comparison-of-speed">&#182;</a></h2><p>we are neither passing, nor returning anything in c++. Just doing same operations in C++ for some array to compare speed.</p>
<p>Mind the overhead for creating new process - it is essential for small sizes.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[13]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">%%</span><span class="k">writefile</span> test_speed.cpp
#include &lt;iostream&gt;
#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;math.h&gt;

// using namespace std;


int main(int argc, char** argv) {
    if (argc &lt; 4){
        std::cout &lt;&lt; &quot;run with n_samples, mean, sigma!&quot;;
        return 1;
    }
    int size = atoi(argv[1]);
    double mean = atof(argv[2]);
    double sigma = atof(argv[3]);
    
    double * data = new double[size];
    double factor = 1. / size;
    for (int i=0; i&lt;size; ++i){
        data[i] = i * factor;
    }
    double result = 0.;
    double s = 0.;
    double x = 0.;
    
    for (int i=0; i&lt;size; ++i){
        x = (data[i] - mean);
        s =  x * x / (2 * (sigma * sigma));
        result += log(exp(-s) / (sqrt(2 * M_PI) * sigma));
    }
    
    std::cout &lt;&lt; std::endl &lt;&lt; result &lt;&lt; std::endl;
    return 0;
}
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>Overwriting test_speed.cpp
</pre>
</div>
</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[14]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">!</span>g++ test_speed.cpp -o test_speed -O3
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[15]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="k">def</span> <span class="nf">llh_cpp</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">mean</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="n">size</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">data</span><span class="p">)</span>
    <span class="n">out</span> <span class="o">=</span> <span class="o">!</span>./test_speed <span class="o">{</span>len<span class="o">(</span>data<span class="o">)}</span> <span class="o">{</span>mean<span class="o">}</span> <span class="o">{</span>sigma<span class="o">}</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Data-generation,-checking-that-all-functions-output-the-same-value">Data generation, checking that all functions output the same value<a class="anchor-link" href="#Data-generation,-checking-that-all-functions-output-the-same-value">&#182;</a></h3>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[16]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">from</span> <span class="nn">collections</span> <span class="kn">import</span> <span class="n">OrderedDict</span>
<span class="n">functions</span> <span class="o">=</span> <span class="n">OrderedDict</span><span class="p">()</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;scipy&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_scipy</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;numpy&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_numpy</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;cython&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_cython</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;numexpr&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_numexpr</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;parakeet&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_parakeet</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;theano&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_theano</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;fortran&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_fortran</span>
<span class="n">functions</span><span class="p">[</span><span class="s">&#39;c++&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">llh_cpp</span>
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[17]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">data</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="mi">1000000</span><span class="p">)</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="s">&#39;float64&#39;</span><span class="p">)</span>
<span class="k">for</span> <span class="n">name</span><span class="p">,</span> <span class="n">function</span> <span class="ow">in</span> <span class="n">functions</span><span class="o">.</span><span class="n">items</span><span class="p">():</span>
    <span class="k">print</span> <span class="n">name</span><span class="p">,</span> <span class="n">function</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="mf">1.1</span><span class="p">)</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>scipy -1431841.8928
numpy -1431841.8928
cython -1431841.8928
numexpr -1431841.8928
parakeet -1431841.8928
theano -1431841.8928
fortran -1431841.90671
c++ None
</pre>
</div>
</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[18]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">timeit</span> 
<span class="n">sizes</span> <span class="o">=</span> <span class="p">[</span><span class="mi">10</span> <span class="o">**</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">10</span> <span class="o">**</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">10</span> <span class="o">**</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">10</span> <span class="o">**</span> <span class="mi">8</span><span class="p">]</span>
<span class="kn">import</span> <span class="nn">pandas</span>
<span class="n">scores</span> <span class="o">=</span> <span class="n">pandas</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">data</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">columns</span><span class="o">=</span><span class="n">functions</span><span class="o">.</span><span class="n">keys</span><span class="p">(),</span> <span class="n">index</span><span class="o">=</span><span class="n">sizes</span><span class="p">)</span>
<span class="k">for</span> <span class="n">size</span> <span class="ow">in</span> <span class="n">sizes</span><span class="p">:</span>
    <span class="k">for</span> <span class="n">name</span><span class="p">,</span> <span class="n">function</span> <span class="ow">in</span> <span class="n">functions</span><span class="o">.</span><span class="n">items</span><span class="p">():</span>
        <span class="n">data</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">size</span><span class="p">)</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="s">&#39;float64&#39;</span><span class="p">)</span>
        <span class="n">result</span> <span class="o">=</span> <span class="o">%</span><span class="k">timeit</span> -o function(data, 0.1, 1.1)
        <span class="n">scores</span><span class="o">.</span><span class="n">loc</span><span class="p">[</span><span class="n">size</span><span class="p">,</span> <span class="n">name</span><span class="p">]</span> <span class="o">=</span> <span class="n">result</span><span class="o">.</span><span class="n">best</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>100 loops, best of 3: 11.2 ms per loop
100 loops, best of 3: 6 ms per loop
100 loops, best of 3: 10.8 ms per loop
100 loops, best of 3: 7.14 ms per loop
100 loops, best of 3: 6.54 ms per loop
100 loops, best of 3: 8.07 ms per loop
100 loops, best of 3: 6.03 ms per loop
100 loops, best of 3: 16.6 ms per loop
10 loops, best of 3: 134 ms per loop
10 loops, best of 3: 62.2 ms per loop
10 loops, best of 3: 108 ms per loop
10 loops, best of 3: 71 ms per loop
10 loops, best of 3: 63.4 ms per loop
10 loops, best of 3: 81.7 ms per loop
10 loops, best of 3: 57.1 ms per loop
10 loops, best of 3: 83 ms per loop
1 loops, best of 3: 1.82 s per loop
1 loops, best of 3: 857 ms per loop
1 loops, best of 3: 1.08 s per loop
1 loops, best of 3: 713 ms per loop
1 loops, best of 3: 660 ms per loop
1 loops, best of 3: 877 ms per loop
1 loops, best of 3: 639 ms per loop
1 loops, best of 3: 694 ms per loop
1 loops, best of 3: 19.8 s per loop
1 loops, best of 3: 9.01 s per loop
1 loops, best of 3: 11.1 s per loop
1 loops, best of 3: 7.12 s per loop
1 loops, best of 3: 6.44 s per loop
1 loops, best of 3: 8.54 s per loop
1 loops, best of 3: 6.39 s per loop
1 loops, best of 3: 6.88 s per loop
</pre>
</div>
</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Results-(time-in-seconds,-less-is-better)">Results (time in seconds, less is better)<a class="anchor-link" href="#Results-(time-in-seconds,-less-is-better)">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[19]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">scores</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[19]:</div>

<div class="output_html rendered_html output_subarea output_execute_result">
<div style="max-height:1000px;max-width:1500px;overflow:auto;">
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>scipy</th>
      <th>numpy</th>
      <th>cython</th>
      <th>numexpr</th>
      <th>parakeet</th>
      <th>theano</th>
      <th>fortran</th>
      <th>c++</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>100000   </th>
      <td>  0.011209</td>
      <td> 0.006000</td>
      <td>  0.010787</td>
      <td> 0.007137</td>
      <td> 0.006539</td>
      <td> 0.008067</td>
      <td> 0.006031</td>
      <td> 0.016596</td>
    </tr>
    <tr>
      <th>1000000  </th>
      <td>  0.134079</td>
      <td> 0.062242</td>
      <td>  0.108327</td>
      <td> 0.071049</td>
      <td> 0.063404</td>
      <td> 0.081651</td>
      <td> 0.057135</td>
      <td> 0.083036</td>
    </tr>
    <tr>
      <th>10000000 </th>
      <td>  1.823849</td>
      <td> 0.857088</td>
      <td>  1.078784</td>
      <td> 0.713240</td>
      <td> 0.659757</td>
      <td> 0.876672</td>
      <td> 0.638656</td>
      <td> 0.694401</td>
    </tr>
    <tr>
      <th>100000000</th>
      <td> 19.826547</td>
      <td> 9.006766</td>
      <td> 11.110369</td>
      <td> 7.124668</td>
      <td> 6.436616</td>
      <td> 8.542589</td>
      <td> 6.389464</td>
      <td> 6.884255</td>
    </tr>
  </tbody>
</table>
</div>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Comparison-to-numpy-time-(less-is-better)">Comparison to numpy time (less is better)<a class="anchor-link" href="#Comparison-to-numpy-time-(less-is-better)">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[20]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">normalized_scores</span> <span class="o">=</span> <span class="n">scores</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
<span class="k">for</span> <span class="n">column</span> <span class="ow">in</span> <span class="n">normalized_scores</span><span class="o">.</span><span class="n">columns</span><span class="p">:</span>
    <span class="n">normalized_scores</span><span class="p">[</span><span class="n">column</span><span class="p">]</span> <span class="o">/=</span> <span class="n">scores</span><span class="p">[</span><span class="s">&#39;numpy&#39;</span><span class="p">]</span>    
</pre></div>

</div>
</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[21]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">normalized_scores</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[21]:</div>

<div class="output_html rendered_html output_subarea output_execute_result">
<div style="max-height:1000px;max-width:1500px;overflow:auto;">
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>scipy</th>
      <th>numpy</th>
      <th>cython</th>
      <th>numexpr</th>
      <th>parakeet</th>
      <th>theano</th>
      <th>fortran</th>
      <th>c++</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>100000   </th>
      <td> 1.868104</td>
      <td> 1</td>
      <td> 1.797846</td>
      <td> 1.189456</td>
      <td> 1.089899</td>
      <td> 1.344443</td>
      <td> 1.005108</td>
      <td> 2.765994</td>
    </tr>
    <tr>
      <th>1000000  </th>
      <td> 2.154161</td>
      <td> 1</td>
      <td> 1.740413</td>
      <td> 1.141503</td>
      <td> 1.018674</td>
      <td> 1.311835</td>
      <td> 0.917955</td>
      <td> 1.334087</td>
    </tr>
    <tr>
      <th>10000000 </th>
      <td> 2.127960</td>
      <td> 1</td>
      <td> 1.258662</td>
      <td> 0.832166</td>
      <td> 0.769766</td>
      <td> 1.022849</td>
      <td> 0.745146</td>
      <td> 0.810186</td>
    </tr>
    <tr>
      <th>100000000</th>
      <td> 2.201295</td>
      <td> 1</td>
      <td> 1.233558</td>
      <td> 0.791035</td>
      <td> 0.714642</td>
      <td> 0.948464</td>
      <td> 0.709407</td>
      <td> 0.764343</td>
    </tr>
  </tbody>
</table>
</div>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Conclusion">Conclusion<a class="anchor-link" href="#Conclusion">&#182;</a></h2><p>Many libraries claim that they can speed up number crunching in python.
Results of this test are floating (+- 0.1), but what we can see</p>
<ol>
<li>numpy turned out to be fastest at moderate sizes of arrays</li>
<li>numpy implementation at least not more complex than others</li>
<li>parakeet was the only to get sensible speed up</li>
</ol>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="Technical-info">Technical info<a class="anchor-link" href="#Technical-info">&#182;</a></h1>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[22]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">multiprocessing</span>
<span class="k">print</span> <span class="n">multiprocessing</span><span class="o">.</span><span class="n">cpu_count</span><span class="p">(),</span> <span class="s">&#39;xeon cores&#39;</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>16 xeon cores
</pre>
</div>
</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[23]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">numpy</span><span class="o">.</span><span class="n">__version__</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[23]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>&apos;1.10.1&apos;</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[24]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">scipy</span><span class="o">.</span><span class="n">__version__</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[24]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>&apos;0.14.0&apos;</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[25]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="kn">import</span> <span class="nn">cython</span>
<span class="n">cython</span><span class="o">.</span><span class="n">__version__</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[25]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>&apos;0.21.1&apos;</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[26]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">numexpr</span><span class="o">.</span><span class="n">__version__</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[26]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>&apos;2.4&apos;</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[27]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">parakeet</span><span class="o">.</span><span class="n">__version__</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[27]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>&apos;0.23.2&apos;</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[28]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="n">theano</span><span class="o">.</span><span class="n">__version__</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt output_prompt">Out[28]:</div>


<div class="output_text output_subarea output_execute_result">
<pre>&apos;0.7.0&apos;</pre>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[29]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">!</span>g++ -v
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>Using built-in specs.
COLLECT_GCC=g++
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/4.6/lto-wrapper
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion=&apos;Ubuntu/Linaro 4.6.3-1ubuntu5&apos; --with-bugurl=file:///usr/share/doc/gcc-4.6/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --program-suffix=-4.6 --enable-shared --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.6 --libdir=/usr/lib --enable-nls --with-sysroot=/ --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --enable-gnu-unique-object --enable-plugin --enable-objc-gc --disable-werror --with-arch-32=i686 --with-tune=generic --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 4.6.3 (Ubuntu/Linaro 4.6.3-1ubuntu5) 
</pre>
</div>
</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[30]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-ipython2"><pre><span class="o">!</span>gfortran -v
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/4.6/lto-wrapper
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion=&apos;Ubuntu/Linaro 4.6.3-1ubuntu5&apos; --with-bugurl=file:///usr/share/doc/gcc-4.6/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --program-suffix=-4.6 --enable-shared --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.6 --libdir=/usr/lib --enable-nls --with-sysroot=/ --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --enable-gnu-unique-object --enable-plugin --enable-objc-gc --disable-werror --with-arch-32=i686 --with-tune=generic --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 4.6.3 (Ubuntu/Linaro 4.6.3-1ubuntu5) 
</pre>
</div>
</div>

</div>
</div>

</div>

<p>
    This post was written in IPython. You can download the notebook from
     <a href='https://github.com/arogozhnikov/arogozhnikov.github.io/tree/master/notebooks'>
    repository</a>.
</p>

